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Abstract 

An efficient linear solver plays an important role while solving partial differential 
equations (PDEs) and partial integro-differential equations (PIDEs) type mathematical 
models. In most cases, the efficiency depends on the stability and accuracy of the 
numerical scheme considered. In this article we consider a PIDE that arises in option 
pricing theory (financial problems) as well as in various scientific modeling and deal 
with two different topics. In the first part of the article, we study several iterative 
techniques (preconditioned) for the PIDE model. A wavelet basis and a Fourier sine 
basis have been used to design various preconditioners to improve the convergence 
criteria of iterative solvers. We implement a multigrid (MG) iterative method. In fact, 
we approximate the problem using a finite difference scheme, then implement a few 
preconditioned Krylov subspace methods as well as a MG method to speed up the 
computation. Then, in the second part in this study, we analyze the stability and the 
accuracy of two different one step schemes to approximate the model. 

Keywords: convolutional integral; preconditioner; stability; convergence. 



1 Introduction 

The pricing of options is a central problem in financial investment. It is important in 
both theoretical and practical point of view since the use of options thrives in the financial 

*The author would like to thank Chris C. Stolk, the KdV institute for Mathematics, University of Ams- 
terdam for introducing him wavelet and Fourier sine preconditioners for elliptic operators. 
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market. In option pricing theory, the study of the Black-Scholes equation is very important 
and interesting (study of a parabolic partial differential equation (PDE)). In recent days, 
researchers have extended the model by looking at the nonlocal effects, which is a linear 
partial integro-differential equation (PIDE). 

We consider such a partial integro-differential equation [?J [9] 

i) " UJ) cau.n. (i) 



dt 



where 



d 2 u(x,t) du(x,t) . . f . 
Cu(x, t) = a — h yu — - ru(x, t) + X J(x - y) (u(y, t) - u(x, t)) dy, 

with initial condition 

u(x, 0) = ip(x), — oo < x < oo. 

Here a > 0, A > 0, with (a, A) ^ (0,0), r > and /i e M, J is the kernel of the model 
and u = u(x,t) represents the option price (contingent claim). A normalized kernel function 
J(x), i.e., J n J(x)dx = 1 has been considered in most of the models [TUl H21 [TT] with suitable 
parameter values. In general, J(x—y) is a kernel function that models the interaction between 
options at positions x and y. The effect of close neighbours x and y is usually greater than 
that from more distant ones; this is incorporated in J. For simplicity we assume that J 
is a non-negative function that satisfies smoothness, symmetry and decay conditions. One 
may consider any J to implement the schemes we discuss in this study. |e~' x ' and y^e -1 ^ 2 
are two sample kernel functions. Boundary conditions are always an issue in these types of 
models. Here one may easily consider BCs [9] 

d 2 u 

— — - = (J, as x —7- ±oo. 
ox 2 

Operator defined by ([T|) with er = 0, // = 0, r = comes while modeling phase transi- 
tions [10], dynamics of neurons in the brain model [HI [H], and population dynamics mod- 
els [17] as well. 

Numerical approximation and analysis of PDEs and PIDEs using finite difference, finite 
element method and the pseudo-spectral method are of ongoing research interest. Specially 
for PIDEs, fast and efficient numerical tools are still to be developed. A clear introduction 
about option pricing models and some finite difference schemes to approximate the models 
can be found in [91 [12] . 

A noble study about the model problem ([TJ) can be found in [T2]. The authors consider a 
European and an American vanilla and barrier options based on the variance gamma process. 
They discuss derivation of (JT|) in detail and approximate the model problem numerically by 
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implementing a finite difference algorithm. They present some numerical experiments on the 
option pricing. But no efficient linear algebra solvers for the discrete equivalent of the model 
as well as the stability and the accuracy analysis of the approximation are discussed. 

In [10], Dugald et. al. consider a nonlocal model of phase transitions of type (1T1) (a = 
0, /i = 0, r = 0). Stability of stationary solution and coarsening of solutions have been 
discussed by the authors. They present a finite element scheme to solve the problem and 
discuss some experimental results. 

In [2] , the author also considers the nonlocal model of phase transitions. He approximates 
the problem using the forward Euler scheme and examines the convergence rate of the scheme. 

A convolutional model of 8 Neuron network has been considered in [3]. The author 
approximates the problem using finite element method in space, then he applies implicit 
schemes for time stepping. Then the author analyzes the error in such an approximation. 

The PIDE model ([1]) is well studied in [7J. They discuss viscosity solution of the model 
followed by a few finite difference approximations. They show that the infinite domain can 
be truncated to a finite domain [A, B] where A and B depend on the decay of the kernel 
function J(x). Thus the problem can be considered as a IBVPs. Considering the kernel of 
the convolution integral as 



and B = —A, where e > is considered so that Jg(x) > e. One may consider the model in 
a spatial periodic domain [2] as well. We use these concepts to approximate the model in a 
finite as well periodic spatial domain. 

There are many other articles those discuss these type of models, but to the best of our 
knowledge the discussion about efficient linear solvers for this type of models is absent. So we 
focus on some fast and efficient numerical schemes as well as the stability and the accuracy 
analysis of two finite difference schemes for the operator acting on ([lj. We start the study 
by approximating the problem using the backward Euler in time for (jTJ) and investigate some 
linear algebra tools to speed up the computational process in Q C E. Then we analyze the 
stability and the accuracy of two different schemes considering Q — M.. 

The article is organized in the following way. We propose and implement several efficient 
linear system solvers to compute solutions of (JT|) in Section [2j Then we discuss the stability 
of an explicit and an semi-implicit scheme in Section El We use Fourier transforms of the 
integro-differential equation for our analysis throughout this study. The accuracy analysis of 
two full discrete schemes as well as a semi-discrete approximation are presented in Section [4] 




the authors in [7] formulate 
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and Section [51 respectively. We finish the article with discussion, conclusions and open 
problems in Section [61 



2 Numerical approximation 

Several standard ordinary differential equation solvers are available and can be used to 
approximate the time derivative. So our main goal, in this study, is to approximate the 
model ([T]) in space domain. Here we first perform a time integration, then look for some fast 
and efficient space integration tools. 

Now one may start with the forward Euler scheme for time stepping (an explicit scheme), 
which uses the values of only previous time step to calculate those of the next. The Algorithm 
is very simple, in that each unknown, at time step n + 1, is calculated independently, so it 
does not require simultaneous solution of equations, and can even be performed easily. But it 
is unstable for large time steps. We have analyzed the stability condition, and the accuracy 
of such a scheme in Section [3] and thereafter. Instabilities are big problems in numerical 
approximation. We want to use large time steps and so we are interested in using implicit 
schemes. 

2.1 An implicit scheme 

We start with the implicit Euler scheme for time integration. We approximate the model 
([[}) in time by 

-Ata d2 ^ n( f"> _ At /2 ^ffi^ + ( 1 + r At) u n (x) - At X [ J(x-y) (u n (y) - u n (x)) dy = u n -\x), 
ox ox Jq 

where u n (x) = u(x,t n ), n > 0. We will demonstrate several schemes to approximate the 
semi-discrete spatial model. For simplicity we write 



where, 



and 



Cu n (x) = A(w n (x)) + C 2 (u n (x)) = u n -\x), (2) 



C^ix)) = -Ata^P- - Atfj,^^- + (1 + rAt)u n (x) 
ox z ox 



C 2 (u n (x)) = -AtX / J(x - y) {u n {y) - u n {x)) dy. 



n 



It is easy to verify that the operator £ acting on ([2]) is an elliptic partial differential opera- 
tor [11]. 
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Now it is our aim to design and implement some fast and efficient solvers for ([2]). We 
start by approximating 

d 2 u n (x) _ U? + i - 2U? + du n {x) _ U? +1 - 



dx 2 h 2 dx 2h 



and 



oo „ 

C 1 u n (x i )= V / J{x l -y){U n {x i )~U n 



(y))dy 



j=-oo 
jV/2-1 



« hJixi-x^^ixi) -U n { Xj )). 

j=-N/2 

Based on these approximations we write the full discrete model as 

Au n = U n-l^ (3) 

which is a system of linear equations with unknown U n . The symbol of the discrete equivalent 
of £ can be written as (see Section |3]) 

A syb {At, h£) =At(l- q(0 + ^ sin 2 - ^ sm(h$ 

Considering 



A syb (At, h£) 

the unknown can be expressed in the Fourier domain as (see Section [3] for details) 

U n (0=9 n (At,hOU°(0- 

Since, 

\A syb (At,hO\ > 1 

for any choice of At and h, the scheme is unconditionally stable (a few discrete symbols of 
this type of operators have been evaluated in detail in next section). 

Now the main difficulty of solving linear systems like is that the maximal eigenvalue 
grows exponentially whereas the minimal eigenvalue is bounded. This situation results in an 
exponential growth of the condition number 

cond(A) = 0(N 2 ) = 0(2 2k ), for some k > 1. 

As a result, any iterative solver becomes slower, and a preconditioning is highly needed. To 
be precise for the Krylov subspace type methods, the solution of the linear system Au = b 
with some uq is 

lk ~^ lu - 2 ( ^y+i ) l|n - n ° lu ' iwu = * Ta *. 



where p(A) is the spectral condition number of A. The convergence of the above expression 
is neat, but it has rarely been presented the convergence of conjugate gradient type methods 
unless p(A) ~ 1 [BJ page 128]. Thus it becomes clear that one needs to find a matrix D such 
that 

B = D-^AD' 1 ' 2 

is well conditioned. It is very popular to replace p(A) in the iterative solvers by p(B), which 
is called preconditioning [HI [IB] and is used for the preconditioned linear system solvers. 
Thus we get the motivation to develop and to compare a few preconditioned solvers based 
on the established and popular preconditioning techniques for local second order elliptic 
operators. We implement and demonstrate the power of multigrid, wavelet as well as Fourier 
preconditioners. Our goal here is to implement preconditioners in a traditional way so that 

1. D is a symmetric and positive definite matrix. 

2. p(B) = 0(1), as N ->■ oo. 

To be specific, we discuss several types of preconditioners below. 

Wavelet Diagonal Preconditioning: One of the most successful preconditioners for el- 
liptic PDEs is the wavelet diagonal preconditioning (WDP) which has been studied in 
details in [181 122], and many other references. Since £ is of elliptic type we attempt 
to implement wavelet diagonal preconditioning to solve (J3J). Suppose £ is defined over 
a periodic domain. Then a preconditioner can be defined by combining two separate 
steps: 

1. Define a basis transformation J 7 (wavelet decomposition operator), given by a 
wavelet transformation, and a wavelet reconstruction operator J 7 * whose columns 
are the elements of the wavelet basis denoted by 

2. Define an invertible diagonal scaling matrix S, whose elements are of the form 
s\ « 2~ 2 ' A ', where |A| denotes the scale index of the wavelet. 

We consider the symmetric proconditioner 

D = F^S 1 ' 2 ? 

as a scaled operator [T8| 122] . Then we define the preconditioned operator (equivalent 
to £ ) by 

rs l/2 Fcrs 1/2 F. 

A detailed discussion about designing such a preconditioner can be found in [18]. A pre- 
conditioner of this type is sensitive with boundary conditions. One may also consider 
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D = ^SJ 7 to define a left or a right preconditioner to implement a preconditioned 
BICG solver. The implementation detail is same as the symmetric preconditioner 
discussed above. 

Fourier Sine Preconditioning: Localization in the position-wave number space is an im- 
portant concept in PDEs and can be extended to PIDEs. Most recently, a frame of 
functions, called windowed Fourier frames, has been employed to solve a variable co- 
efficient second order elliptic PDE [3]. Here it is our aim to design and to implement 
preconditioners based on the Fourier sine transformation (FSP) for the PIDE This 
preconditioning is sensitive with boundaries, and works very well for periodic boundary 
value problems. 

The symbol of the operator C defined by can be written as 

Ata£ 2 - AtfiiC + (1 + rAt) - AtAv / 2^(J(0 - J(0)). (4) 

When £ is very large, £ 2 term becomes the dominating term in (J4]) and so C can be 
approximated by Ata^ 2 in the frequency domain. Thus we approximate 

d 2 

Cu « Ata^-^u = Ata^lb k sm(£ k x) . 



dx' 



Let M k = Ata(i ^ 0, then 



Thus 



Now 



i d 2 



l . d 2 

——Atu ——zU 

M k dx 2 



y^bn sin(£ fe x). 



^2 h sin(£ fc x) 



^ h sin(£ 



kX, 



u 



where B is the frame upper bound [16]. It is clear from [16] that B < oo if we consider 
a tight frame. Thus 



i d 2 

——Ata—^u 

M k dx 2 



< oo, and so 



< OO. 



Based on this idea we define a preconditioned operator 
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with the symmetric preconditioner 

P = F*M~ l/2 F, 

where F stands for the Fourier sine transformation operator |16j . The invertibility of 
the operator P£P, considering F a windowed Fourier transform operator, has been 
proved in [I] (where £ is a second order elliptic operator). The idea can be extended 
to the PIDE model that we have considered here. So we avoid attempting to prove the 
invertibility of the Fourier sine preconditioner (FSP) PCP here. One may also consider 
P = F*M~ 1 F to define a left or a right preconditioner to implement a preconditioned 
BIGG solver. The implementation detail is same as the symmetric precomnditioner we 
have discussed above. 

Multigrid Preconditioning: Multigrid (MG) methods are now a days the fastest and 
most efficient numerical solvers for linear systems. There are huge recent literatures 
on MG methods. Actually multigrid method combines two separate ideas [6| [21] : 

1. fine grid residual smoothing by relaxation. 

2. coarse grid residual correction. 

Here the idea is to perform a few iterations (smoothing) in a fine grid, then switch 
to a coarser level and perform a few iterations, and so on. This is called coarse grid 
corrections. After corrections, one switches back to the fine grid and performs a few 
post-smoothing. Thus a multigrid algorithm uses three basic and old steps: 

• relaxation step. 

• restriction step. 

• interpolation step. 

A detailed discussion about multigrid can be found in [51 El [2TJ and in many other 
references. Since the operator (j2J) is of elliptic type, multigrid would be one of the 
choices to be considered to verify it's efficiency. Here we implement a so-called v— cycle 
to solve the system fl3}. It behaves well with both periodic and non-periodic boundary 
conditions. 

In our problem we use just one f-cycle. One can use v w-cycles if the solution is not 
sufficiently accurate after the completion of one cycle. We follow the Algorithm [1] for com- 
putation. 

Algorithm 1. Multigrid method to solve system of linear equations 

To solve the system of linear equations Au = f using the multigrid method 
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INPUT the finest grid matrix A h , right-hand vector f h , L the number of steps to travel 
down the coarsest grid, \x the number of relaxation (iterations) on each grid, tolerance 
Tol, number of v-cycles v, initial solution u = 

OUTPUT The approximate solutions u h . 

Step 1 For 7 = 1, 2, ■ ■ ■ ,v or error < Tol do the following steps 

Step 2 Relax A h u h = f h /x times using the Jacobi iteration with the initial data Mq. 

Step 3 Set r h = f h - A h u h . 



Step 4 For k = 2, 3, • • • , L - 1 



define residual f kh = r_ kh , where r\ 



Mi — S k ~ l ) h 



N 



it 



kh 



0, and relax A hh u kh = f kh fi times as in step 2. 



take the initial guess 



Step 5 Set f h = r Lh and solve A Nh u Lh = f^ h exactly. 
Step 6 For k = L — 1, L — 2, • • ■ ,1 we upgrade u kh by using 



u 



kh 



kh 
2j-l 

kh 
U 2 j 



<_x + «f +1)ft ; j = i,2,---,iv (fc+1)h , 



<+2 



(k+l)h (k+l)h 



J = 1,2, 



, N( k+ i) h - 1, 



ii--, 



it 



kh 
2N ( , 



+ 



U 



(k+l)h (k+l)h 



+ U 



N, 



(k+l)h 



and upgrade the solution \x times as in step 2. 

Step 7 If 1 1?/ 4 1 1 < Tol or for v > 7 > 1 if \\u hr/ — u h,,y+1 \\ < Tol, output the required solution 
u h else "program stopped after v v- cycle" . 



STOP 



This is to note that one may use FFT for each matrix vector multiplication to reduce 
the computational costs since the operator A acting on (JHJ) is a toeplitz matrix [T3] . 



2.2 Numerical results and discussions 

Here we present some experimental/computer generated results to demonstrate the efficiency 
of the schemes. We implement the schemes in MATLAB. The MATLAB function "FFT" is 
used to define the Fourier sine preconditioner; MATLAB functions "wavedec" , and " waverec" 
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have been used for the wavelet diagonal preconditioner with the Daubechies wavelet 'db6'. 
Here we consider a spatial periodic [0 1] domain and 



J{x) 



oo 

E 

r=— oo 



J°°(x-r) with J°°(x) 



100 



7T 



A detailed discussion about such a consideration of the kernel function can be found in [3] . 
We consider o = 0.01, p = 0.01, r = 0.01, A = 0.1, 5 = 100 for all the numerical results 
presented here. 

In Figure [U we present condition numbers of the preconditioned operators PCP, and 
DCD, as well as the condition number of C. We notice that p(DCD), and p(PCP) are of 
(9(1), where as p(C) is of 0(N 2 ). Then, in Figure [2], we compare the number of iterations 
taken by the preconditioned solvers for a set of N values. We notice that the preconditioned 
systems converge in a few iterations and the number of iterations is independent of the 
system size. Then we demonstrate the total CPU time taken to solve the linear system 



5.5 



4.5 



3.5 



* * * * 



-fc-wavelet-precond. 
□ FS-precond. 



100 200 300 400 500 600 700 800 900 1000 
N 



10" 



10 



10 



10 



- -0 



100 200 300 400 500 600 700 800 900 1000 
N 



Figure 1: Left Figure: Condition numbers of the wavelet preconditioned operator, and the 
Fourier sine preconditioned operator, both are of 0(1), Right Figure: Condition number of 
A, which is of 0(N 2 ). 



by the solvers MG, WDP CG and FSP CG respectively to see the time efficiency of the 
techniques in Figure [3j Here we observe that in terms of CPU time the MG out performs all 
other schemes. In fact, the MG method takes very little computational time compared to the 
other two. The WDP and FSP techniques take most of the time to define the preconditioners, 
the preprocessing steps to use preconditioned linear system solvers. 



10 




Figure 2: The Number of iterations taken to converge by the conjugate gradient, the WDP 
conjugate gradient and the FSP conjugate gradient methods to solve (j2J) considering At = 
0.01. 



2.3 An explicit implicit scheme 

While solving the linear system ([3]) we notice that A is a full matrix. Thus matrix vector 
multiplications are computationally costly. To reduce the computation cost further we look 
for an another scheme that may reduce computational costs. We implement an explicit 
implicit scheme where A becomes a spare matrix, thus reduces computation costs in matrix 
vector multiplications. 

We approximate the model ([1]) in time by 

-Ata \ ' +{l+rAt)u n {x) = u n -\x)+Atfi +At\ / J{x-y) {u n ~\y) - u n -\x)) dy, 

ox ox Jq 

where u n (x) = u(x,t n ), n>0. For simplicity we write 

C l {u n {x)) = C 2 u n -\x), (5) 



where 

d{u n {x)) = -Ata - I } + (1 + rAt)u n {x), 

ox z 

and 

dv n ~ 1 (r) r 
£ 2 (u n (x)) = Atfi y ' + AtX / J(x - y) (u n_1 (y) - u n ~\x)) dy. 
ox J n 
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Figure 3: CPU time taken to converge by the WDP, the FSP and the MG methods for 
various choices of system size. For the multigrid method we consider one log 2 (TV) — 2 level 
v— cycle with one SOR iteration with u — 1.2 in all levels but the coarsest one. In the 
coarsest level we solve the system exactly. We use Intel(R) Core(TM) i3 CPU M380 at 2.53 
GHz processor with ram 2.00 GB. Here we solve ([3]) considering At = 0.01, and by varying 
spatial grid points N, and uq(x) = exp[— 100(x — .5) 2 ]. 



The operator L\ is an elliptic partial differential operator [TTJ. After the time integration, 
the right hand side of (J5J) is a known vector and explicitly depends on u"" 1 , n > 1. Thus all 
the linear algebra tools we discussed above for (j2J) are applicable to (jSJ), and they are indeed, 
efficient schemes for elliptic PDEs. 

To justify our claim we implement the MG method, the fastest tool we implemented in 
the previous section, to solve the linear system obtained from (0). We compare the CPU 
time taken to solve the linear system obtained by the implicit solver and the explicit implicit 
solver (J3J in Figure HI Here we notice that the scheme ([2]) and the explicit implicit scheme 
(JSJ) are comparable. In fact, it is observed from Figure H] that the scheme (jSJ) requires a 
minimum CPU time to converge compared to all other solvers. 



3 Stability analysis 

From the above Section we see that the scheme (jSJ) dominates the implicit scheme in terms 
of computational time. This numerical experiment motivates us to analyze the stability and 



12 




Figure 4: CPU time taken to converge by MG methods for various choices of system size 
from the implicit as well as the explicit implicit solvers. Here we consider one log 2 (TV) — 2 
level v— cycle with one SOR iteration with u = 1.2 in all levels but the coarsest one. In 
the coarsest level we solve the system exactly. We use Intel(R) Core(TM) i3 CPU M380 
at 2.53 GHz processor with ram 2.00 GB. Here we solve (|3]) and ([5]) considering At = 0.01, 
and by varying spatial grid points N, and uq(x) = exp[— 100(x — .5) 2 ]. 



the accuracy of an explicit and an explicit implicit scheme. For the simplicity of the stability 
analysis we consider a — \i — \ — r — 1. Here we analyze the stability of the forward Euler 
scheme (explicit) and a mix Euler scheme (explicit implicit). We consider the linear partial 
integro-differential equation [9] (an IVP) 

/oo 
J(x - y) (u(y,t) - u(x,t))dy (6) 
-oo 

with u(x,to) = Uo(x), i6R. This IVP can be approximated in space by 

dUjjt) U j+1 - 2Uj + gj_i U j+1 - Uj-i 
dt h 2 2h j 

oo 

+ h J{xj-x k ){U k -Uj) (7) 

k=— oo 

for each j G Z where Uj(t) ~ u(xj,t) and Xj = jh where h is the uniform spacing between 
the grid points Xj and Xj+i for all j G Z. We need the following definitions to support our 
study. 
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For the sequence {v m : m G Z} on the mesh points {x m = mh : m G Z} the discrete 
Fourier Transform (DFT) is defined by 



2lX 



e- ihmi v r , 



m=—oo 



if w m G L 2 (/iZ), and its inverse is 



where £ G [-^, ?]. Parseval's Formulae [191 [23] are defined as 

oo 

i*koi 2 #= e %»r = 



~ 1 1 2 



(9) 



(10) 



An explicit scheme 

We apply the explicit Euler scheme to the semi-discrete model ([7]) to obtain 

Jjn _ 7V™ r/n _ OTjn i rm 



jjn+i _ u: 



J 



2h 



h 2 



+ hAt Jfa ~ x k ) {U r k l - U^) 



k=— oo 



where C/j 1 = U(xj,t n ). This is equivalent to 

7V™ _ Tjn Jjn _ njjn , jjn 



At- 



2h ' W 

oo 

+hAt J2 J ( x j ~ x k) u k- 

k=—oo 



f/j 1 ( 1 - At - hAt J ( x j 



- x k 



k=—oo 



'IF 



We multiply ( TlT]) by -7=e and sum over all j to obtain 



, oo , oo 

E e~ ijh( UV+ 1 = e ~ iMu j [l- At -hAt J ( x i - 



^7T 



J=-oo 



+ 



2^ 



j=-oo 
oo 



fc=— oo 



2^ 

h 



/27T 



=— oo 

OO / 

e-^(At 



j=-oo 
oo 



hAt J ( x 3 ~ x k )e- i(j - kM 

j=-oo 



/> 2 



OO / 

At 



2vr . 

j=-oo 



Jjn _ Jjn 

2h 
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So using J(x) = J(—x) we have 



Thus 



where 



U n+1 (0 = I 1 - At + hAt J2 J ( x 3 - x k) {e t{k ~ m - 1) \ U n (0 

[_ j=-oo J 



U n (£) = (9(h^At)TU (Oi (12) 



g(h£, At) = 1 - At + At I h ^ e^J(x r ) - h e~ irm J{ 

\ r=— oo r=— oo 

+ ^ (e** + e-^ - 2) + ^ (e«* - e"**) 
= l-At + v^FAt(j(0- J(0)) -4^sin 2 ^ + ^sin(^). (13) 

Now we carry out the stability analysis of fill I) following [TJ [TJ] . We need the following 
Lemma to bound g(h£, At). 

Proposition 1. Assume that J(x) G L 2 (M) D C(R) satisfies 
HI J(x) > 0; 

H2 J(x) zs normalized such that J(x)dx = 1; 
H3 J(x) zs symmetric, i.e. J(x) = J(—x), for all x G R; 
H4 J(x) zs decreasing on (0, oo); 
H5 J(0 > 0. 



TTjen HI - H4 give the DFT results < J(0) and J(0 < J(0) < Jf + J(0 /or a// 

£ G [-f,f] and t/ie CFT resn/ts J(f) < J(0) < + J(£). Further, if H5 ftoZtfe, t/ien 
J(f) > /or a// J G # r (M), r > ±, Jf. 

Now we back to the main discussion. The scheme is stable if 

\g\ < i- 
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Here 



h 2 V 2 



\g\ 2 = (l- At + V2^At(J(0- J(0))-^sin 2 ' ^ 
+ (^Ysin 2 (/>0<1 



■"> 



gives 



A f ( 1 + ^)-*) + l^(^)-?f^(f, + 
16.^4 A«\ , 1 ■ 2,, A „ A, ,8. 2 ^ 



M Sin ifJ+^ Sin ^JM 2 - 2 « (f)+ /? Sin 
Thus applying Proposition [TJ we have 

. / 25 16\ 

and so 

t - (3/i 2 + 4) 2 + /i 2 - (3/i 2 + 4) 2 ~ (3 + 4//i 2 ) 2 ' ( > 



Theorem 1. If J(x) is a normalized symmetric nonnegative function and J G L 2 (K)nC(I 

trik < \\U°\\ h 



then there exists < ( 3+4 / fe 2)2 < At* stzc/i that 



for all < At < At* and n>0. 

Proof. The proof easily follows from perse val's relation. □ 

Thus in the discrete L2 norm, ( Till is a stable scheme [191 Definition 1.5.1] with the 
stability condition (JTj 



An explicit implicit scheme 

Applying a mixed Euler scheme we write a full discrete version of the model (jUJ) by 



C/n+^fjn = _ AtU n+i + M ^+i L + A t^ ^zl 

00 

k=—oo 

1(3 



where £7™ = U(xj,t n ). This is equivalent to 

Tjn _ Tjn TJ n+l — 9TT n 

f/" +1 (l + At) = At j+ \ U * +At i+1 7, 
J 2h h z 

oo 

+hAt J ( x i- x k)U%. 



( 1 - hAt J(xj - x k ) 

(15) 



fc=— oo 



fc=— oo 



Multiplying (j!5p by -J=e ij7 ^ and summing over all j we get 



(i 



J=-OC 



< ' jhi l'j [l-hAt J ( x j - x k) 



2tt . 

j=-oo 



h 
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h 



J2 ' im u% 



k=— oo 
oo 



fc=— oo 
oo 



ZiAt J(xj — Xfc)e" 



j=-oo 



2tt . 

j=-oo 



v 7 ^ ^ V h J 

v j=— oo x ' 



So using J(x) = J(—x) 
At 



f/ n+1 (0(l + At - — + - 2 )) = { 1 + hAt J2 J ( X J ~ x k) (e i{k ~ m - l) }> U n (0 



j=-oo 



giving 

And we write 
where 

g(h£ t At) 
The scheme is stable if 



+ (e<* - 1) 



u n+1 (0 = g(h^At)u n (0- 

U n (0 = (g(hZ,At)TU°(0, 



1 + V2^At - J(0)) + f (e ih S - 1) 



1 + At + 4f sin 2 f 



(16) 



(17) 



1 + v^At (j(£) - J(0)) + ^ (e th « - 1) 



< 



At . 2 h£ 
1 + At + 4— sin 2 

nr 2 
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which gives 

Now 
and 



2ttA* ( J(0 - J(0) 



At 
~h 



At 
~h 



< 



At 



At 



4—- sin 

h 2 2 



At 

h 2 



27T J(0" J(0) 



At) 2 < 1 + At 



< 



sin 



2 

< 2. 



Simplifying the above inequality we get 

2 2 



At 2 <f(0 - 1 



g(0 < At 2 



2 



29(0 



and so 



At < 



2h(h 



3h 2 + 4h + 2' 

Theorem 2. If J(x) is a normalized symmetric nonnegative function and J G L 2 { 
then there exists < 3^2+4^2 — ^* suc ^ ^at 

\\U n \\h < \\U°\\h 

for all < At < At* and n>0. 

Proof. The proof easily follows from perse val's relation. 



nc(i 



□ 



Thus in the discrete L 2 norm, (I15|) is a stable scheme [T9J Definition 1.5.1] with the 
stability condition (jl8p . We demonstrate maximum values of At from both ( TT4"j) and 018p 
respectively in Figure |3] for various choices of /i. It shows the dominance of the semi- implicit 
scheme. 



Computational algorithm 

From the schemes f[T2l and fTTol) it follows that the DFT gives a each way to compute 
numerical solutions. The approximate solution U n (-) can be computed in the spatial domain 
simply, accurately and rapidly using the following steps. For faster computations, one may 
precompute the FFT of u , J(x), and J(0). 

1. Compute the fast Fourier transform (FFT) of Uq. 

2. Compute g using the FFT of J. 

3. Evaluate g n and multiply with the result in step 1. 

4. Compute the inverse FFT of the product defined in step 3. 



18 



0.2 



0.4 0.6 
h 



0.8 



Figure 5: Maximum choices of At from the inequalities 014p and (I18j) . 



4 Accuracy analysis 

Applying the continuous Fourier transform ([6]) can be written as 

ut(U) = q(0HU), (19) 

where 

VV Z7T \/ ZlX \/ ZlX / 

Thus the exact solution of (1191) in the frequency domain is 

w(£,t)=e^ (0- (20) 

Here it is easy to verify that 9ft(g) < (Proposition (TJ which is presented in Section |3]) which 
gives the stability property |tt(£,t)| < |{to(£, 

Computational algorithm 

The following steps can be taken to compute the exact solution and the error in schemes 
(H2D and (US]). 

1. Compute the FFT of uq. 
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2. Compute q as denned in (|T9|) using FFT of J. 



3. Evaluate exp(nAtq) and multiply with the result obtained from step 1. 

4. Compute the inverse FFT of the product defined in 3. 

5. Evaluate \\u(-,t) - U n (-)\\. 

In this section it is our aim to present a theoretical bound of the error term \\u — U n \\. 
Now we carry out the convergence analysis of (ITT]) and (fl5l) following [19]. We apply the 
inverse CFT on (120]) to get 



which is the exact solution of (jH]). 

4.1 The explicit scheme (1121) 

Using the inverse DFT formula OH]) on (fl2"|) . the approximate solution can be presented as 



Applying the Fourier interpolation (19] the mesh function ( 122]) can be written as 




(21) 




(22) 




(23) 



Thus 



u(x,t n )-SU n (x) 





(24) 



So 



u(x,t n )-SU n (x)\\ 2 < 




e^u (0-(9(hZ^t)) n MO\ 2 dZ 




(25) 



using Parseval's relation and the stability property q < 0. 
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Let us find a bound related to the-evolution error first. Here 



< 



2n 
2 



lfl<£ 



l£l<i 



V W|£|< f ^ V W 



since At) | < 1. Now following pEl page 204], [2] 

27TJ 



l£l<i 



5>({ 



^ < d^^iitxoii^ 



(26) 



where C\(a) = 2 (\) 2 ° Y^?=i (2j — 1) 2 ° assuming that the initial function is smooth and 



there exists a > ~ such that ||mo||_h" ct (m) is bounded, and 

V 27T 7|f|>? 



(27) 



'ICl>i 

When t n = nAt 



n-l 



r=0 



Since < and \g(h£, At)| < 1 we have 



or equivalently 

J c 

Now, for the scheme ffTTT) 



\e^ Atn - g n \ < n\e^ At 



g(ht,At) n \<n\e^ At -g(h(,At)\. 



(28) 



<?(/*£, At) = e 



, «f 



2tt V2ir 



1 - At + At\Fhx ( J(f ) - J(0) 



At . 2 ^ zAt . 
— 4— sin — H — — sin{n^) 



h 2 



h 



AW2-K J(f ) - J(0) - 



e 2 + i ^ 



27T V27T 



AtV2vr J(f ) - J(0) 



4 . 2 />£ , iAt . 
_= — sm — —stn(hf) — 1 



J'=2 



+E-r h^F i(£)-i(o) 



e 2 + i + ^ 



^7T V27T 



(29) 
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Assuming that J G H r (§L) with r > ~ and applying the Poisson summation formula, (|29|) 
becomes 

= -A t ^Wi« + ^)-^)) + ^ S in=(f)-A^ 

j^O ^ J \ J 

iAt 

+iAt£ - — sin (h£) + 0(At 2 ) 



h 



• At 

' (h0 3 + O((h£f) + O(At 2 ). (30) 



Proposition 2. |1J/ Assume that HI, H3 and H5 of Lemma U\ hold and in addition, the 
following condition holds: 

H6. |J(O<0/or^>0. 

Then, for all \£\ < = |E jy o + ¥) " | < 2 ^(f )■ 

Thus applying Proposition [21 (|30|) can be written as 

\e Am) - g(h£,At)\ < Atd(h) + C 2 Ath 2 \£\ A (31) 



where d(h) = 2V2^/(f). If J G L 2 (M), then |J(f)| -» as |f | -> oo [H], [201 page 30]. 
The rate of convergence determines the accuracy of the scheme. We have 



f \(e^-g(h{,Aty)u (0\ 2 dt 
J \£\<% 

< [ n 2 \e^ At -g(ht,At)\ 2 \u (O\ 2 d^ using (EH]) 

[ \Atd(h) + C 2 Ath 2 \£\ 3 \ 2 \u (0\ 2 d^ using (JSU) 

/oo 
|Ci(/i) + c a /i 2 |e| 4 | 2 M0l a de 
-oo 



< n 2 



< 

J — oo 

< tnC^WuoW 2 + t n C 2 h 2 \\u \\ 2 H2m . (32) 
Thus applying lj2Ej) . (l27j) and (l3"2]h (120] takes the form 

\\u(x,t n )-SU n (x)\\ < C 1 (h)\\u \\ 2 + C 2 h\\u \\l 2m +C 3 (o-)h°\\u \\ H « m (33) 
for all Mo G H a (M.) with a > ~. Thus we end up with the following result. 
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Theorem 3. If the kernel function J(x) satisfies assumptions HI - H6 and < f77]j is a stable 
approximation for the IDE then there exist constants C\(h), C2, Cs(a) such that 

\\u(x,Q -SU n (x)\\ < t n C x {h)\\u \\ +C 2 h\\u \\ H , m + C 3 (a)h°\\uo\\H°(R) 
for any uq G H a (M) with a > |. 



4.2 The explicit implicit scheme (1151) 

Using series expansion 



e \ ^ = l + Atv 7 ^ I J(0 - J(0) - + 

V v2vr v27T/ 



+ At 2 ( V2n ( J(0 - J(0) - ^ + -|=) ) 2 + C(At 



Also 

<?(X, At) = (l + v^At (J(0 - J(0)) + ^ (e«« - 1) Vl + At + 4^ sin 2 ^ * . 
Letting At < — . 1 „ n = — — fc2 , he < 1 (since minesin 2 ^ = 0), considering 4ri and ^ 

to — l+Asin 2 ^- /i 2 +4sin 2 ^ — v P y ' to h ' ft 2 

constants we have 

, . At . 2 /iA _1 , /. At . 2 /iA /. At . 2 hC 7 
1 + At + 4— sin 2 — = 1 - At + 4— sin 2 — + At + 4— sin 2 



/i 2 2 J \ h 2 2 J \ h 2 2 

and so 

, A * :„2^ _1 \+zf£\ ( 1 , , A t . 



Atg(0 x ^1 + At + 4— sin 2 -yl = Atg(£) ^1 - \At + 4— sin 2 ) 
where q(g) = ^ - J(0)V Also 

^ (e** - 1) = ^ " ^ 2 e 2 /2 - ^) + 0{h^% 



gives 



At ( i/i 3 £ 3 \ / ( . At h£ 
1 ih£ - h 2 i 2 j2 2_ 1 _ At + 4— sin 2 s 



h V eyvV ^ 2 2 



,2,2 /9 ^ 3 e 3 \ 4 At 2 . 2 ^ 

-(^-^/2- — U— sin y . 
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Thus 



g(h£,At) 



1 + At 



1 - T7T 




'f + m + 1 Ui - hr/2 



+At 2 



h 



1 



( 



ih£ - h 2 £ 2 /2 





+ 0(At 3 ) 



gives 



e 



Atv^l J(£)-J(0) 




g(h£, At) < Atd(h) + C 2 Ath 2 ^ + C 3 At 2 . 



Thus following similar procedure as of the accuracy analysis of the explicit Euler scheme we 
estimate the accuracy of the scheme (I15p by the following theorem. 

Theorem 4. If the kernel function J(x) satisfies assumptions HI - H6 and (T73]) is a stable 
approximation for the IDE (0|) ; then there exist constants Ci(h), C 2 , C%, C^a) such that 

\\u(x,t n ) - SU n (x)\\ < t n Ci(h)\\u \\ + C 2 h\\u \\ H 2 m + C 3 At||u || + C 4 (cr) /i ct ||m ||^(k), 

for any u G H a (M) with & > \. 

We compute error in such approximations that have been presented above considering 
various choices of the kernel function and the initial function. We present errors estimated 
by Theorem [3J and Theorem H] in Figure El From this computation we observe the supremacy 
of the explicit implicit scheme as well as the importance of the choices of the initial function 
uq(x) and the kernel function J{x). Here it can easily be noticed that smooth J(x) and 
u (x) give better accuracy and that justifies the Theorem [3] and the Theorem HI 

5 Accuracy of the semidiscrete approximation 

Here we study the accuracy of the scheme (J7j) . Applying the discrete Fourier transform on 



(34) 



where £ G [-- , \ ] and 
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t=0.5, ||u(:,t)-U n (:)|| 



t=0.5, ||u(:,t)-U n (:)|| 




Figure 6: Here we present the error ||w(:, t n ) — U n (:)\\ estimated by the explicit scheme (right 
two) and the explicit-implicit (left two). Here in the bottom Figures we consider J(x) = e~' x ', 
uq(x) = e - '^; in the top Figures we consider Uq(x) = J— e~ Wx \ J(x) = 



and thus 

Applying the inverse Fourier transform to (135]) 

7T 

1 ft .-, 



U m (t) 

We interpolate U m (t) defined in fl36|) bv [T9] 

1 

SU{x,t) 



2vr J-e 

h 



(35) 



(36) 
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Similar to the Theorem [3] (using (12T]) ). 
u(x, t) — SU(x, t) - 



2,71 



(0 - eW (0 ) dC 



Itt 



(37) 



Theorem 5. If J, and J satisfy the assumptions HI - H6 and Uq G H a (R) with a > |, 
t/ien i/iere exzsi constants Ci(h), C^{a) such that 

\\u{x,t) -SU(t)\\ < tCi(h)\\uo\\ + C2/i||«o||fl2(R) + C 3 (a)h a \\u \\ H ^(^), 
where |7]) is a semidiscrete approximation to the IDE |7^). 
Proof. We have 



\u{;t)-SU(;t)\\' 



\u{-,t) -SU(-,t)\ 2 dx 



< 



1 



Itt 
1 



iei<f 



^7T 



(38) 



lel>f 



since by Lemma [T]i?eaZ(g) < 0. Similar to the analysis of the full discrete approximation the 
first part of the right-hand side of (|38p can be written as 



1 



< 



/2tt J\{\<« 
2 



«£ (£)-e«W (£) 



di 



'2u J\£\< 



e^ t -e^ t f\u (t)\ 2 d£ + 



2tt J\z\< 



X>(« 



27TJ 



We have 



| e «(0t_ e 9(€)t| <t|g(0-g(OI 
since rea/(g(£)) < and rea/(g(£)) < 0. Now 



9(f) - 5(f) 



27TJ 



2^ ^ 



h 
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- J 

- J 



z sin(/i£) 

27TJ 

h 



Thus 



and C(h) = 2J(f ) as h -> 0. Now 

f \ e mt _ e mf \u (C)\ 2 dC < e [ m) - g(0| 2 \u (0\ 2 dC 

J \t\<i J \z\<l 

<t 2 C(h) 2 \\u \\ 2 + C 2 h 2 \\u \\ 2 H2m , 

gives 

I \e mt - e*<«f |fio(e)| a de < t 2 C 2 (h)\\u f + C 2 t 2 h 2 \\u \\ 2 Hm . (39) 
J \t\<% 

Thus applying fl2"H|l .(l2"7jl and (155), (I3"g|) takes the form 

||m(z,*) - SU(t)\\ < tCj^llMoll + C 2 h\\u \\ 2 Hm + C 3 (a)h a \\u \\ H « m (40) 
for some d, C 2 , C 3 for all u e H a (R) with a > \. □ 

6 Summary and conclusions 

In this study, we consider a linear partial integro-differential operator (PIDO) that comes in 
modeling financial engineering problems as well as in modeling various scientific problems. 
We study a few finite difference schemes (FDSs) for European style options with a jump- 
diffusion term (the PIDO). In the first part of the study we introduce several preconditioned 
linear system solvers for the full discrete equivalent of the model. We observe that all 
the preconditioned solvers are very efficient, and the multigrid solver is way better than the 
wavelet diagonal preconditioned solver and the Fourier sine preconditioned solvers. In fact, a 
one v— cycled Multigrid solver is several times faster than the other two. The implementation 
costs for the sine and the wavelet preconditioning are relatively higher than that of the 
multigrid technique. So we conclude that a multigrid method can be used to speed up the 
computation of the finite dimensional (full discrete) PIDE model. Here we also conclude 
that the explicit implicit scheme outperforms the implicit scheme in terms of computational 
costs. 

Here, in the second part of this study, we analyze the stability and the accuracy of two 
different finite difference schemes. While analyzing the stability and the accuracy of the 
finite difference schemes (an explicit scheme as well as an explicit implicit scheme) we notice 
that the schemes are conditionally stable (under some reasonable restrictions imposed on 
the kernel function). The explicit implicit scheme is faster than that of the explicit scheme 
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as well as the implicit scheme, which agrees with the properties of the time and the space 
discretizations of the PIDE we consider in this study. We establish some bounds of the error 
in such full discrete as well as semi-discrete schemes. 

Here we analyze the model in one space dimension only. Preconditioners can be employed 
to speed up the computational process for the full discrete model, specially for two and three 
space dimensional domains as well as preconditioned solvers along with higher order multi- 
step schemes may be better options to think of, and that leaves as future research directions. 
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